###################################################################################################
# Appellate Court Influence over District Courts in the United States #############################
# Michael Olson and Albert Rivero #################################################################
# Analysis ########################################################################################
###################################################################################################

###################################################################################################
# working directory and packages ##################################################################
###################################################################################################

  library(lfe)
  library(stargazer)
  library(tidyverse)
  library(ggplot2)
  library(readxl)
  library(viridis)
  library(interflex)
  library(fixest)
  library(english)
  library(stringr)
  library(DataCombine)

  setwd("C:/Users/michael.p.olson/Dropbox/Research/Appellate court response to Supreme Court/full_replication_package")
  
###################################################################################################
# use panel means data to demonstrate variability of the measure and how it changes over time #####
################################################################################################### 
  
# for all circuits over time
  
  read_rds("panel_means.rds") %>%
    filter(month == 1,
           circuit != "Federal") %>%
    mutate(mean_panel = -1 * mean_panel,
           median_all = -1 * median_all,
           circuit = factor(circuit, levels = c("First",
                                                "Second",
                                                "Third",
                                                "Fourth",
                                                "Fifth",
                                                "Sixth",
                                                "Seventh",
                                                "Eighth",
                                                "Ninth",
                                                "Tenth",
                                                "Eleventh",
                                                "D.C."))) %>%
    select(`Mean Panel Liberalism` = mean_panel, `Circuit Median` = median_all, year, circuit) %>%
    pivot_longer(-c(year, circuit)) %>%
    rename(Measure = name,
           Value = value,
           Year = year) %>%
    ggplot(aes(x = Year, y = Value, linetype = Measure)) +
    geom_line() +
    facet_wrap(~ circuit) +
    scale_x_continuous(breaks = c(1970, 1990, 2010)) +
    theme_classic()+
    theme(legend.position = "bottom",plot.margin = unit(c(1, 0.25, 0.25, 0.25), "cm"))
  
  ggsave("../new_measure_comparison.jpeg",width=8,height=5.2,units="in", dpi = 800)
  
# for the DC circuit in the 1980s  
  
  pm <- read_rds("panel_means.rds") %>%  filter(circuit == "D.C.",year == 1985 | year==1986 | year==1987 | year==1984)
  pm <- pm %>%  mutate(mean_panel = -1 * mean_panel,
                       median_all = -1 * median_all) %>%
    select(`Mean Panel Liberalism` = mean_panel, `Circuit Median` = median_all, year, circuit, month)
  pm <- pm %>% pivot_longer(-c(year, month, circuit)) %>%
    rename(Measure = name,
           Value = value,
           Year = year,
           Month = month)
  pm$Month <- zoo::as.yearmon(paste(pm$Year, pm$Month), "%Y %m")
  
  jcs_circuit_month <- read_csv("jcs_circuit_month.csv")
  
  s <- jcs_circuit_month[jcs_circuit_month$circuit=="D.C." & (is.element(jcs_circuit_month$yeara,c(1984:1987)) | (is.element(jcs_circuit_month$yearl,c(1984:1987)))),] %>% arrange(yeara,mona,yearl,monl)
  s[s$monl==12,"yearl"] <- s[s$monl==12,"yearl"]+1
  s$monl <- s$monl+1
  s[s$monl==13,"monl"] <- 1
  s$yearc <- s$yeara
  s$yearc[s$yeara<1984] <- s$yearl[s$yeara<1984]
  s$monc <- s$mona
  s$monc[s$yeara<1984] <- s$monl[s$yeara<1984]
  s$modec <- "Arrives"
  s$modec[s$yeara<1984] <- "Departs"
  s$ln <- unlist(lapply(str_split(s$judge,","),function(x) x[1]))
  s$event <- paste(s$ln,s$modec)
  s <- s[,c("yearc","monc","event")]
  s$Month <- zoo::as.yearmon(paste(s$yearc, s$monc), "%Y %m")
  
  s <- left_join(s,pm[pm$Measure=="Mean Panel Liberalism",])
  s <- s[s$Month=="Oct 1985" | s$Month=="Oct 1986",]
  s <- s %>% group_by(Month) %>% summarize(event=paste(event,collapse=",\n",sep=", "))
  s$y <- c(-0.15,-0.5)
  s$Month2 <- s$Month
  s$Month[1] <- s$Month[1]-1
  s$Month[2] <- s$Month[2]+0.75
  
  ggplot(data=pm,aes(x = Month, y = Value, linetype = Measure)) +
    geom_line(size=1.5) +
    xlab("Year")+
    geom_text(data=s,aes(x=Month,y=y,label=event),inherit.aes = F)+
    geom_segment(aes(x=s$Month[1],y=s$y[1]+0.075,xend=s$Month2[1],yend=pm$Value[pm$Measure=="Circuit Median" & pm$Month==s$Month2[1]]-0.01),
                 arrow = arrow(length = unit(0.5, "cm")),show.legend = F)+
    geom_segment(aes(x=s$Month[1],y=s$y[1]+0.075,xend=s$Month2[1]-0.03,yend=pm$Value[pm$Measure=="Mean Panel Liberalism" & pm$Month==s$Month2[1]]-0.01),
                 arrow = arrow(length = unit(0.5, "cm")),show.legend = F)+
    geom_segment(aes(x=s$Month[2]-0.35,y=s$y[2],xend=s$Month2[2]+.1,yend=pm$Value[pm$Measure=="Circuit Median" & pm$Month==s$Month2[2]]+0.04),
                 arrow = arrow(length = unit(0.5, "cm")),show.legend = F)+
    geom_segment(aes(x=s$Month[2]-0.35,y=s$y[2],xend=s$Month2[2]+0.04,yend=pm$Value[pm$Measure=="Mean Panel Liberalism" & pm$Month==s$Month2[2]]-0.01),
                 arrow = arrow(length = unit(0.5, "cm")),show.legend = F)+
    scale_colour_grey(start=0,end=0.6)+
    theme_classic()+
    theme(legend.position = "bottom",plot.margin = unit(c(1, 0.25, 0.25, 0.25), "cm"))
  
  ggsave("../dc_example.jpeg",width=8,height=5.2,units="in", dpi = 800)
  
###################################################################################################
# Create example of the 10th Circuit ##############################################################
###################################################################################################
  
  jcs_circuit_month <- read_csv("jcs_circuit_month.csv") %>%
    filter(circuit == "Tenth") %>%
    mutate(JCS = -1 * JCS,
           JCS = JCS + 0.4767033,
           JCS = JCS / 1.056703)
  
  # We will grab all possible panels per month and take the median
  
  my_tibble <- tibble(year = rep(1965:2017, each = 12),
                      month = rep(1:12, times = length(1965:2017)))
  panel_combn <- list()
  jcs <- list()
  
  for(i in 1:636) {
    i_year <- my_tibble$year[i]
    i_month <- my_tibble$month[i]
    
    # We'll grab a) all judges who arrive STRICTLY before i_year and depart after
    # i_year, then b) for judges who arrive in i_year, those who arrive on/after
    # the current month, and c) for judges who depart in i_year, those who depart
    # on/before the current month
    
    i_jcs_a <- jcs_circuit_month %>%
      filter(yeara < i_year, yearl > i_year) %>%
      select(JCS, judge)
    
    i_jcs_b <- jcs_circuit_month %>%
      filter(yeara == i_year, mona <= i_month, yeara!=yearl) %>%
      select(JCS, judge)
    
    i_jcs_c <- jcs_circuit_month %>%
      filter(yearl == i_year, monl >= i_month, yeara!=yearl) %>%
      select(JCS, judge)
    
    i_jcs_d <- jcs_circuit_month %>%
      filter(yeara == i_year, yearl == i_year, mona <= i_month, monl >= i_month) %>%
      select(JCS, judge)
    
    i_jcs <- bind_rows(i_jcs_a, i_jcs_b, i_jcs_c, i_jcs_d)
    
    if(length(i_jcs %>% pull(JCS)) > 3) {
      panel_combn[[i]] <- combn(i_jcs %>% pull(JCS), 3) %>% apply(2, median)
    } else if(length(i_jcs) %in% c(1:3)) {
      panel_combn[[i]] <- i_jcs %>% pull(JCS) %>% median()
    }
    jcs[[i]] <- i_jcs %>% pull(JCS)
  }
  
  my_tibble <- my_tibble %>%
    mutate(panel_combn = panel_combn,
           jcs = jcs)
  
  combns <- tibble(JCS = my_tibble %>% 
                     filter(year == 2012, month == 12) %>%
                     pull(panel_combn) %>%
                     unlist()) %>%
    count(JCS)
  
  example <- jcs_circuit_month %>%
    filter(yeara < 2012, yearl > 2012) %>%
    select(JCS, judge)
  
  example2 <- jcs_circuit_month %>%
    filter(yeara == 2012, mona <= 12) %>%
    select(JCS, judge)
  
  example3 <- jcs_circuit_month %>%
    filter(yearl == 2012, monl >= 12) %>%
    select(JCS, judge)
  
  example <- bind_rows(example,
                       example2,
                       example3)
  
  for_plot <- example %>%
    mutate(judge = str_extract(judge, "^([^,]*)")) %>%
    group_by(JCS) %>%
    summarize(Judges = paste(judge, collapse = ", ")) %>%
    mutate(Judges = paste("(", Judges, ") ", round(JCS, 3), sep = "")) %>%
    left_join(combns) %>%
    mutate(n = ifelse(is.na(n), 0, n),
           n = n / sum(n))
  
  panel_mean <- my_tibble %>% 
    filter(year == 2012, month == 12) %>%
    pull(panel_combn) %>%
    unlist() %>%
    mean()
  
  my_breaks <- c(for_plot %>% slice(1:3) %>% pull(JCS),
                 panel_mean,
                 for_plot %>% slice(4:8) %>% pull(JCS))
  
  my_labels <- c(for_plot %>% slice(1:3) %>% pull(Judges),
                 expression(italic("MEAN: 0.303")),
                 for_plot %>% slice(4:8) %>% pull(Judges))
  
  for_plot %>%
    ggplot(aes(x = JCS, y = n, fill = JCS)) + 
    geom_col(show.legend = FALSE) +
    scale_x_continuous(breaks = my_breaks, labels = my_labels) +
    scale_y_continuous(labels = scales::percent_format(accuracy = 1)) +
    scale_fill_gradient(low = "darkred",
                        high = "darkblue") + 
    geom_vline(xintercept = panel_mean,
               size = 1,
               linetype = "dashed") +
    theme_classic() +
    coord_flip() +
    xlab("Liberalism") +
    ylab("Percentage of Panels as Median")
  ggsave("../CA10.jpeg", width = 8, height = 6, units = "in", dpi = 800)
  
###################################################################################################
# load main analysis data #########################################################################
###################################################################################################

  district <- read_csv("district.csv")
  
###################################################################################################
# show variability of our main independent variable ###############################################
###################################################################################################  
  
  circ_year <- aggregate(mean_panel~circuit+year,data=district,FUN=mean)
    circ_year <- dplyr::arrange(circ_year,circuit,year)
    circ_year <- DataCombine::slide(Var="mean_panel",GroupVar="circuit",TimeVar="year",slideBy=-1,data=circ_year)
    circ_year <- DataCombine::slide(Var="mean_panel",GroupVar="circuit",TimeVar="year",slideBy=-3,data=circ_year)
    circ_year <- DataCombine::slide(Var="mean_panel",GroupVar="circuit",TimeVar="year",slideBy=-5,data=circ_year)
    circ_year <- DataCombine::slide(Var="mean_panel",GroupVar="circuit",TimeVar="year",slideBy=-10,data=circ_year)
    circ_year$change1 <- abs(circ_year$mean_panel-circ_year$`mean_panel-1`)
    circ_year$change3 <- abs(circ_year$mean_panel-circ_year$`mean_panel-3`)
    circ_year$change5 <- abs(circ_year$mean_panel-circ_year$`mean_panel-5`)
    circ_year$change10 <- abs(circ_year$mean_panel-circ_year$`mean_panel-10`)
    
    circ_year <- circ_year[,grep("mean\\_panel",colnames(circ_year),invert=T)]
    circ_year <- circ_year %>% gather(Years,change,change1:change10)
    circ_year$Years <- str_to_title(as.english(as.numeric(gsub("change","",fixed=T,circ_year$Years))))
    circ_year$Years <- factor(circ_year$Years,levels=c("One","Three","Five","Ten"))
    
    jpeg("../year_change.jpeg",width=8,height=5, units = "in", res = 800)
    print(  
      ggplot(circ_year,aes(x=change,group=Years,colour=Years))+
        stat_ecdf(size=2)+
        scale_colour_viridis_d(end=0.9)+
        theme_minimal()+
        ylab("P(Absolute Change < X)")+
        xlab("Absolute Change In 'Mean Panel Liberalism' Within Circuit Courts")
    )
    dev.off()

  jpeg("../year_change_density.jpeg",width=8,height=5, units = "in", res = 800)
  print(  
    ggplot(circ_year,aes(x=change,group=Years,colour=Years,fill=Years))+
      geom_density(alpha=0.3)+
      scale_colour_viridis_d(end=0.9)+
      scale_fill_viridis_d(end=0.9)+
      theme_minimal()+
      ylab("Density")+
      xlab("Absolute Change In 'Mean Panel Liberalism' Within Circuit Courts")
  )
  dev.off()
    
###################################################################################################
# district court responsiveness to the mean liberalism of panels ##################################
###################################################################################################

# Correlation between circuit liberalism and district judge voting

  panel_1 <- felm(libcon ~ mean_panel | 0 | 0 | circuit, data = district)

# Let's try within judge

  panel_2 <- felm(libcon ~ mean_panel | judge_fe | 0 | circuit, data = district)

# Adding year fixed effects

  panel_3 <- felm(libcon ~ mean_panel | judge_fe + factor(year) | 0 | circuit, data = district)

# Adding casetype fixed effects

  panel_4 <- felm(libcon ~ mean_panel | judge_fe + factor(year) + factor(casetype) | 0 | circuit,
                  data = district,
                  exactDOF = TRUE)

# Judge-category fixed effects

  panel_5 <- felm(libcon ~ mean_panel | judge_category + factor(year) + factor(casetype) | 0 | circuit,
                  data = district,
                  exactDOF = TRUE)

# make table

  stargazer(panel_1, panel_2, panel_3, panel_4, panel_5,
            style = "AJPS",
            keep.stat = c("n", "adj.rsq"),
            header = FALSE,
            covariate.labels = c("Mean Liberalism of Circuit Panels"),
            label = "felm",
            title = "Predicting liberal district court voting with the mean liberalism of appellate panels; coefficients from linear fixed effect models (standard errors clustered by circuit)",
            digits = 2,
            model.names = FALSE,
            add.lines = list(c("Judge Fixed Effects", "", "\\checkmark", "\\checkmark", "\\checkmark", ""),
                             c("Judge-Category Fixed Effects", "", "", "", "", "\\checkmark"),
                             c("Year Fixed Effects", "", "", "\\checkmark", "\\checkmark", "\\checkmark"),
                             c("Case Type Fixed Effects", "", "", "", "\\checkmark", "\\checkmark")),
            dep.var.labels = "District Court Judge Liberal Vote") %>%
    cat(sep = "\n", file = "../felm.tex")

  
  # Add Democratic Presidential Vote Share
  
  panel_1_dem <- felm(libcon ~ mean_panel + dem_share | 0 | 0 | circuit, data = district)
  
  # Let's try within judge
  
  panel_2_dem <- felm(libcon ~ mean_panel + dem_share  | judge_fe | 0 | circuit, data = district)
  
  # Adding year fixed effects
  
  panel_3_dem <- felm(libcon ~ mean_panel + dem_share  | judge_fe + factor(year) | 0 | circuit, data = district)
  
  # Adding casetype fixed effects
  
  panel_4_dem <- felm(libcon ~ mean_panel + dem_share  | judge_fe + factor(year) + factor(casetype) | 0 | circuit,
                  data = district,
                  exactDOF = TRUE)
  
  # Judge-category fixed effects
  
  panel_5_dem <- felm(libcon ~ mean_panel + dem_share  | judge_category + factor(year) + factor(casetype) | 0 | circuit,
                  data = district,
                  exactDOF = TRUE)
  
  
  # make table
  
  stargazer(panel_1_dem, panel_2_dem, panel_3_dem, panel_4_dem, panel_5_dem,
            style = "AJPS",
            keep.stat = c("n", "adj.rsq"),
            header = FALSE,
            covariate.labels = c("Mean Liberalism of Panels","Dem. Pres. Vote Share"),
            label = "felm_dem",
            title = "Predicting liberal district court voting with the mean liberalism of appellate panels, controlling for Democratic presidential vote share; coefficients from linear fixed effect models (standard errors clustered by circuit)",
            digits = 2,
            model.names = FALSE,
            add.lines = list(c("Judge Fixed Effects", "", "\\checkmark", "\\checkmark", "\\checkmark", ""),
                             c("Judge-Category Fixed Effects", "", "", "", "", "\\checkmark"),
                             c("Year Fixed Effects", "", "", "\\checkmark", "\\checkmark", "\\checkmark"),
                             c("Case Type Fixed Effects", "", "", "", "\\checkmark", "\\checkmark")),
            dep.var.labels = "Liberal Vote") %>%
    cat(sep = "\n", file = "../felm_dem.tex")
  
  # Test senior judge measure
  
  panel_1_sen <- felm(libcon ~ mean_panel_senior | 0 | 0 | circuit, data = district)
  panel_2_sen <- felm(libcon ~ mean_panel_senior | judge_fe | 0 | circuit, data = district)
  panel_3_sen <- felm(libcon ~ mean_panel_senior | judge_fe + factor(year) | 0 | circuit, data = district)
  panel_4_sen <- felm(libcon ~ mean_panel_senior  | judge_fe + factor(year) + factor(casetype) | 0 | circuit,
                  data = district,
                  exactDOF = TRUE)
  panel_5_sen <- felm(libcon ~ mean_panel_senior | judge_category + factor(year) + factor(casetype) | 0 | circuit,
                  data = district,
                  exactDOF = TRUE)
  
  stargazer(panel_1_sen, panel_2_sen, panel_3_sen, panel_4_sen, panel_5_sen,
            style = "AJPS",
            keep.stat = c("n", "adj.rsq"),
            header = FALSE,
            covariate.labels = c("Mean Liberalism of Panels"),
            label = "felm_senior",
            title = "Predicting liberal district court voting with the mean liberalism of appellate panels (senior judges included); coefficients from linear fixed effect models (standard errors clustered by circuit)",
            digits = 2,
            model.names = FALSE,
            add.lines = list(c("Judge Fixed Effects", "", "\\checkmark", "\\checkmark", "\\checkmark", ""),
                             c("Judge-Category Fixed Effects", "", "", "", "", "\\checkmark"),
                             c("Year Fixed Effects", "", "", "\\checkmark", "\\checkmark", "\\checkmark"),
                             c("Case Type Fixed Effects", "", "", "", "\\checkmark", "\\checkmark")),
            dep.var.labels = "Liberal Vote") %>%
    cat(sep = "\n", file = "../felm_senior.tex")
  
# different types of standard errors/inference
  
# state-clustered 

  st_panel_1 <- felm(libcon ~ mean_panel | 0 | 0 | state, data = district)
  
  st_panel_2 <- felm(libcon ~ mean_panel | judge_fe | 0 | state, data = district)
  
  st_panel_3 <- felm(libcon ~ mean_panel | judge_fe + factor(year) | 0 | state, data = district)
  
  st_panel_4 <- felm(libcon ~ mean_panel | judge_fe + factor(year) + factor(casetype) | 0 | state,
                  data = district,exactDOF = TRUE)
  
  st_panel_5 <- felm(libcon ~ mean_panel | judge_category + factor(year) + factor(casetype) | 0 | state,
                  data = district,exactDOF = TRUE)
  
# district-clustered

  di_panel_1 <- felm(libcon ~ mean_panel | 0 | 0 | statdist, data = district)
  
  di_panel_2 <- felm(libcon ~ mean_panel | judge_fe | 0 | statdist, data = district)
  
  di_panel_3 <- felm(libcon ~ mean_panel | judge_fe + factor(year) | 0 | statdist, data = district)
  
  di_panel_4 <- felm(libcon ~ mean_panel | judge_fe + factor(year) + factor(casetype) | 0 | statdist,
                     data = district,exactDOF = TRUE)
  
  di_panel_5 <- felm(libcon ~ mean_panel | judge_category + factor(year) + factor(casetype) | 0 | statdist,
                     data = district,exactDOF = TRUE)
  
# load wild bootstrap results from stata

  circuit_wild <- as.matrix(read_excel("circuit_wild_ci.xlsx",col_names=F))
  
  state_wild <- as.matrix(read_excel("state_wild_ci.xlsx",col_names=F)) 
  
  district_wild <- as.matrix(read_excel("district_wild_ci.xlsx",col_names=F)) 
  
# combine into data for plot
  
  coefs <- c(coef(panel_1)[2],coef(panel_2)[1],coef(panel_3)[1],coef(panel_4)[1],coef(panel_5)[1])
  
  circuit_se <- c(sqrt(panel_1$clustervcv[2,2]),sqrt(panel_2$clustervcv[1,1]),sqrt(panel_3$clustervcv[1,1]),
                sqrt(panel_4$clustervcv[1,1]),sqrt(panel_5$clustervcv[1,1]))
  circuit_ci <- cbind(coefs-1.96*circuit_se,coefs+1.96*circuit_se)
  
  state_se <- c(sqrt(st_panel_1$clustervcv[2,2]),sqrt(st_panel_2$clustervcv[1,1]),sqrt(st_panel_3$clustervcv[1,1]),
                sqrt(st_panel_4$clustervcv[1,1]),sqrt(st_panel_5$clustervcv[1,1]))
  state_ci <- cbind(coefs-1.96*state_se,coefs+1.96*state_se)
  
  district_se <- c(sqrt(di_panel_1$clustervcv[2,2]),sqrt(di_panel_2$clustervcv[1,1]),sqrt(di_panel_3$clustervcv[1,1]),
                sqrt(di_panel_4$clustervcv[1,1]),sqrt(di_panel_5$clustervcv[1,1]))
  district_ci <-cbind(coefs-1.96*district_se,coefs+1.96*district_se)

  ci_dat <- as.data.frame(rbind(
                  state_ci,
                  circuit_ci,
                  district_ci,
                  state_wild,
                  circuit_wild,
                  district_wild))

  rownames(ci_dat) <- NULL; colnames(ci_dat) <- c("ci_lower","ci_upper")

  ci_dat$Type <- factor(c(rep("State-Clustered SE",5),rep("Circuit-Clustered SE",5),rep("District-Clustered SE",5),
                   rep("State-Blocked Wild Bootstrap",5),rep("Circuit-Blocked Wild Bootstrap",5),rep("District-Blocked Wild Bootstrap",5)))

  ci_dat$Model <- factor(rep(c(1:5),6))
  ci_dat$coef <- rep(coefs,6)
  
# make plot

  jpeg("../alt_ses.jpeg",width=8,height=5, units = "in", res = 800)
  print(
    ggplot(ci_dat,aes(x=Model,y=coef,ymin=ci_lower,ymax=ci_upper,colour=Type,group=Type))+
      geom_linerange(position=position_dodge(width=0.5),size=2)+
      geom_point(colour="black",size=4)+
      geom_hline(yintercept=0,colour="black",size=2,linetype=2)+
      ylab("Estimate of Effect of 'Mean Panel Liberalism'")+
      xlab("Model (Corresponding to Table 1)")+
      scale_colour_viridis_d()+
      theme_minimal()
  )
  dev.off()
  
###################################################################################################
# rate of appellate review ########################################################################
###################################################################################################
  
  app_rev_1 <- felm(libcon ~ mean_panel * appellate_rev  | 0 | 0 | circuit, 
                  data = district %>% filter(year >= 1971),
                  exactDOF = TRUE)
  
  app_rev_2 <- felm(libcon ~ mean_panel * appellate_rev  | judge_fe  | 0 | circuit, 
                  data = district %>% filter(year >= 1971),
                  exactDOF = TRUE)
  
  app_rev_3 <- felm(libcon ~ mean_panel * appellate_rev  | judge_fe + year | 0 | circuit, 
                  data = district %>% filter(year >= 1971),
                  exactDOF = TRUE)
  
  app_rev_4 <- felm(libcon ~ mean_panel * appellate_rev  | judge_fe + year + casetype | 0 | circuit,
                  data = district %>% filter(year >= 1971),
                  exactDOF = TRUE)
  
  app_rev_5 <- felm(libcon ~ mean_panel * appellate_rev  | judge_category + year + casetype | 0 | circuit, 
                  data = district %>% filter(year >= 1971),
                  exactDOF = TRUE)
  
# write out the coefficient on the interaction from the main model

  fileConn<-file("../int_coef.tex")
  writeLines(format(round(coef(app_rev_4)[4],3)), fileConn)
  close(fileConn)
  
# use the interflex package to make a marginal effects plot

  review_binning <- interflex(data=as.data.frame(district[district$year>=1971,]),
                              estimator = "binning",
                              Y="libcon",D="mean_panel",X="appellate_rev",
                                  FE=c("judge_fe","year","casetype"),cl="circuit",
                                  na.rm=T,nbins = 4,full.moderate = F,
                                  vcov.type = "cluster",
                                  ylab = "Marginal Effect of \'Mean Panel Liberalism\'",
                                  xlab = "Moderator: Rate of Appellate Review",
                                  theme.bw = T)
  
  jpeg("../review_me.jpeg",width=8,height=5, units = "in", res = 800)
  plot(review_binning,
       theme.bw=T,bin.labs=F,
       xlab="Rate of Appellate Court Review",
       ylab="Marginal Effect of \'Mean Panel Liberalism\'")
  dev.off()
  
  
# rate of reversal
  
  app_reverse_1 <- felm(libcon ~ mean_panel * rate_reversal | 0 | 0 | circuit, 
                    data = district %>% filter(year >= 1971 & district$rate_reversal<.1),
                    exactDOF = TRUE)
  
  app_reverse_2 <- felm(libcon ~ mean_panel * rate_reversal | judge_fe  | 0 | circuit, 
                    data = district %>% filter(year >= 1971 & district$rate_reversal<.1),
                    exactDOF = TRUE)
  
  app_reverse_3 <- felm(libcon ~ mean_panel * rate_reversal | judge_fe + year | 0 | circuit, 
                    data = district %>% filter(year >= 1971 & district$rate_reversal<.1),
                    exactDOF = TRUE)
  
  app_reverse_4 <- felm(libcon ~ mean_panel * rate_reversal | judge_fe + year + casetype | 0 | circuit,
                   data = district %>% filter(year >= 1971 & district$rate_reversal<.1),
                   exactDOF = TRUE)
  
  app_reverse_5 <- felm(libcon ~ mean_panel * rate_reversal | judge_category + year + casetype | 0 | circuit, 
                    data = district %>% filter(year >= 1971 & district$rate_reversal<.1),
                    exactDOF = TRUE)
  
  
# reversal figure (excluding E.D.Tex. 1993)

  reversal_binning <- interflex(data=as.data.frame(district[district$year>=1971 & district$rate_reversal<.1,]),
                                estimator = "binning",
                                Y="libcon",D="mean_panel",X="rate_reversal",
                                  FE=c("judge_fe","year","casetype"),cl="circuit",
                                  na.rm=T,nbins = 4,full.moderate = F,
                                  vcov.type = "cluster",
                                  ylab = "Marginal Effect of \'Mean Panel Liberalism\'",
                                  xlab = "Moderator: Rate of Reversal",
                                  theme.bw = T)
  
  jpeg("../reversal_me.jpeg",width=8,height=5, units = "in", res = 800)
  plot(reversal_binning,
       theme.bw=T,bin.labs=F,
       xlab="Rate of Appellate Court Reversal",
       ylab="Marginal Effect of \'Mean Panel Liberalism\'")
  dev.off()
  
# make table for review
  
  stargazer(app_rev_1, app_rev_2, app_rev_3, app_rev_4, app_rev_5,
            style = "AJPS",
            keep.stat = c("n", "adj.rsq"),
            header = FALSE,
            covariate.labels = c("Mean Liberalism of Panels",
                                 "Rate of Appellate Review",
                                 "Mean Liberalism $\\times$ Rate of Review"),
            label = "appellate_rev_tab",
            title = "Predicting liberal district court voting with the mean liberalism of appellate panels, interacted with a measure of the rate of appellate review for the previous calendar year. The data are limited to years after 1970. Coefficients are from linear fixed effect models with standard errors clustered by circuit.",
            digits = 2,
            model.names = FALSE,
            add.lines = list(c("Judge Fixed Effects", "", "\\checkmark", "\\checkmark", "\\checkmark", ""),
                             c("Judge-Category Fixed Effects", "", "", "", "", "\\checkmark"),
                             c("Year Fixed Effects", "", "", "\\checkmark", "\\checkmark", "\\checkmark"),
                             c("Case Type Fixed Effects", "", "", "", "\\checkmark", "\\checkmark")),
            dep.var.labels = "Liberal Vote") %>%
    cat(sep = "\n", file = "../felm_rev.tex")
  
# make table for reversal
  
  stargazer(app_reverse_1, app_reverse_2, app_reverse_3, app_reverse_4, app_reverse_5,
            style = "AJPS",
            keep.stat = c("n", "adj.rsq"),
            header = FALSE,
            covariate.labels = c("Mean Liberalism of Panels",
                                 "Rate of Appellate Reversal",
                                 "Mean Liberalism $\\times$ Rate of Reversal"),
            label = "appellate_reverse_tab",
            title = "Predicting liberal district court voting with the mean liberalism of appellate panels, interacted with a measure of the rate of appellate reversal for the previous calendar year. The data are limited to years after 1970. Coefficients are from linear fixed effect models with standard errors clustered by circuit.",
            digits = 2,
            model.names = FALSE,
            add.lines = list(c("Judge Fixed Effects", "", "\\checkmark", "\\checkmark", "\\checkmark", ""),
                             c("Judge-Category Fixed Effects", "", "", "", "", "\\checkmark"),
                             c("Year Fixed Effects", "", "", "\\checkmark", "\\checkmark", "\\checkmark"),
                             c("Case Type Fixed Effects", "", "", "", "\\checkmark", "\\checkmark")),
            dep.var.labels = "Liberal Vote") %>%
    cat(sep = "\n", file = "../felm_reverse.tex")

  # Rate of Review with judge party

  # Check judge_fe works:
  district %>% count(judge_fe, party) %>% select(-n) %>% count(judge_fe) %>% filter(n> 1)

  app_party_rev_1 <- felm(libcon ~ mean_panel * appellate_rev + I(party==1) * appellate_rev | 0 | 0 | circuit,
                          data = district %>% filter(year >= 1971, multiple_parties == 0),
                          exactDOF = TRUE)

  app_party_rev_2 <- felm(libcon ~ mean_panel * appellate_rev + I(party==1):appellate_rev | judge_fe  | 0 | circuit,
                          data = district %>% filter(year >= 1971, multiple_parties == 0),
                          exactDOF = TRUE)

  app_party_rev_3 <- felm(libcon ~ mean_panel * appellate_rev + I(party==1):appellate_rev | judge_fe + year | 0 | circuit,
                          data = district %>% filter(year >= 1971, multiple_parties == 0),
                          exactDOF = TRUE)

  app_party_rev_4 <- felm(libcon ~ mean_panel * appellate_rev + I(party==1):appellate_rev | judge_fe + year + casetype | 0 | circuit,
                          data = district %>% filter(year >= 1971, multiple_parties == 0),
                          exactDOF = TRUE)

  app_party_rev_5 <- felm(libcon ~ mean_panel * appellate_rev + I(party==1):appellate_rev | judge_category + year + casetype | 0 | circuit,
                          data = district %>% filter(year >= 1971, multiple_parties == 0),
                          exactDOF = TRUE)

  stargazer(app_party_rev_1, app_party_rev_2, app_party_rev_3, app_party_rev_4, app_party_rev_5,
            style = "AJPS",
            keep.stat = c("n", "adj.rsq"),
            header = FALSE,
            covariate.labels = c("Mean Liberalism of Panels",
                                 "Rate of Appellate Review",
                                 "Judge is a Democrat",
                                 "Mean Liberalism $\\times$ Rate of Review",
                                 "Judge is a Democrat $\\times$ Rate of Review"),
            label = "appellate_rev_party_tab",
            title = "Predicting liberal district court voting with the mean liberalism of appellate panels,
            interacted with a measure of the rate of appellate review for the previous calendar year which is also interacted with an indicator for the judge's party.
            The data are limited to years after 1970. Coefficients are from linear fixed effect models with standard errors clustered by circuit.",
            digits = 2,
            model.names = FALSE,
            add.lines = list(c("Judge Fixed Effects", "", "\\checkmark", "\\checkmark", "\\checkmark", ""),
                             c("Judge-Category Fixed Effects", "", "", "", "", "\\checkmark"),
                             c("Year Fixed Effects", "", "", "\\checkmark", "\\checkmark", "\\checkmark"),
                             c("Case Type Fixed Effects", "", "", "", "\\checkmark", "\\checkmark")),
            dep.var.labels = "Liberal Vote") %>%
    cat(sep = "\n", file = "../felm_rev_judgeparty.tex")


  # Rate of Reversal with judge party

  app_party_reverse_1 <- felm(libcon ~ mean_panel * rate_reversal + I(party==1) * rate_reversal | 0 | 0 | circuit,
                              data = district %>% filter(year >= 1971, multiple_parties == 0),
                              exactDOF = TRUE)

  app_party_reverse_2 <- felm(libcon ~ mean_panel * rate_reversal + I(party==1):rate_reversal | judge_fe  | 0 | circuit,
                              data = district %>% filter(year >= 1971, multiple_parties == 0),
                              exactDOF = TRUE)

  app_party_reverse_3 <- felm(libcon ~ mean_panel * rate_reversal + I(party==1):rate_reversal | judge_fe + year | 0 | circuit,
                              data = district %>% filter(year >= 1971, multiple_parties == 0),
                              exactDOF = TRUE)

  app_party_reverse_4 <- felm(libcon ~ mean_panel * rate_reversal + I(party==1):rate_reversal | judge_fe + year + casetype | 0 | circuit,
                              data = district %>% filter(year >= 1971, multiple_parties == 0),
                              exactDOF = TRUE)

  app_party_reverse_5 <- felm(libcon ~ mean_panel * rate_reversal + I(party==1):rate_reversal | judge_category + year + casetype | 0 | circuit,
                              data = district %>% filter(year >= 1971, multiple_parties == 0),
                              exactDOF = TRUE)

  stargazer(app_party_reverse_1, app_party_reverse_2, app_party_reverse_3, app_party_reverse_4, app_party_reverse_5,
            style = "AJPS",
            keep.stat = c("n", "adj.rsq"),
            header = FALSE,
            covariate.labels = c("Mean Liberalism of Panels",
                                 "Rate of Appellate Reversal",
                                 "Judge is a Democrat",
                                 "Mean Liberalism $\\times$ Rate of Reversal",
                                 "Judge is a Democrat $\\times$ Rate of Reversal"),
            label = "appellate_reverse_party_tab",
            title = "Predicting liberal district court voting with the mean liberalism of appellate panels,
            interacted with a measure of the rate of appellate reversal for the previous calendar year which is also interacted with an indicator for the judge's party.
            The data are limited to years after 1970. Coefficients are from linear fixed effect models with standard errors clustered by circuit.",
            digits = 2,
            model.names = FALSE,
            add.lines = list(c("Judge Fixed Effects", "", "\\checkmark", "\\checkmark", "\\checkmark", ""),
                             c("Judge-Category Fixed Effects", "", "", "", "", "\\checkmark"),
                             c("Year Fixed Effects", "", "", "\\checkmark", "\\checkmark", "\\checkmark"),
                             c("Case Type Fixed Effects", "", "", "", "\\checkmark", "\\checkmark")),
            dep.var.labels = "Liberal Vote") %>%
    cat(sep = "\n", file = "../felm_reverse_judgeparty.tex")


###################################################################################################
# heterogeneity by category #######################################################################
###################################################################################################
  
  panel_cat_1 <- felm(libcon ~ mean_panel * factor(category) | 0 | 0 | circuit, data = district)
  
  panel_cat_2 <- felm(libcon ~ mean_panel * factor(category) | judge_fe | 0 | circuit, data = district)
  
  panel_cat_3 <- felm(libcon ~ mean_panel * factor(category) | judge_fe + factor(year) | 0 | circuit, data = district)
  
  panel_cat_4 <- felm(libcon ~ mean_panel + mean_panel : factor(category) | judge_fe + factor(year) + factor(casetype) | 0 | circuit,
                      data = district,
                      exactDOF = TRUE)
  
  panel_cat_5 <- felm(libcon ~ mean_panel + mean_panel : factor(category) | judge_category + factor(year) + factor(casetype) | 0 | circuit,
                      data = district,
                      exactDOF = TRUE)

# make table 
    
  stargazer(panel_cat_1, panel_cat_2, panel_cat_3, panel_cat_4, panel_cat_5,
            style = "AJPS",
            keep.stat = c("n", "adj.rsq"),
            header = FALSE,
            covariate.labels = c("Mean Liberalism of Panels",
                                 "Civil Rights/Civil Liberties",
                                 "Economics",
                                 "Mean Liberalism $\\times$ CR/CL",
                                 "Mean Liberalism $\\times$ Economics"),
            label = "felm_category",
            title = "Predicting liberal district court voting with the mean liberalism of appellate panels, interacted by category of case.  The omitted category is criminal cases.  Coefficients are from linear fixed effect models with standard errors clustered by circuit.",
            digits = 2,
            model.names = FALSE,
            add.lines = list(c("Judge Fixed Effects", "", "\\checkmark", "\\checkmark", "\\checkmark", ""),
                             c("Judge-Category Fixed Effects", "", "", "", "", "\\checkmark"),
                             c("Year Fixed Effects", "", "", "\\checkmark", "\\checkmark", "\\checkmark"),
                             c("Case Type Fixed Effects", "", "", "", "\\checkmark", "\\checkmark")),
            dep.var.labels = "Liberal Vote") %>%
    cat(sep = "\n", file = "../felm_category.tex")
  
###################################################################################################
# heterogeneity by workload #######################################################################
###################################################################################################

  workload_1 <- felm(libcon ~ mean_panel * logwork  | 0 | 0 | circuit, 
                     data = district,
                     exactDOF = TRUE)
  
  workload_2 <- felm(libcon ~ mean_panel * logwork  | judge_fe | 0 | circuit, 
                     data = district,
                     exactDOF = TRUE)
  
  workload_3 <- felm(libcon ~ mean_panel * logwork  | judge_fe + year| 0 | circuit, 
                     data = district,
                     exactDOF = TRUE)
  
  workload_4 <- felm(libcon ~ mean_panel * logwork  | judge_fe + year + casetype | 0 | circuit, 
       data = district,
       exactDOF = TRUE)

  workload_5 <- felm(libcon ~ mean_panel * logwork | judge_category + year + casetype | 0 | circuit, 
                     data = district,
                     exactDOF = TRUE)
  
# make binned figure

  workload_binning <- interflex(data=as.data.frame(district),Y="libcon",D="mean_panel",X="logwork",
                                estimator = "binning",
                                  FE=c("judge_fe","year","casetype"),cl="circuit",
                                  na.rm=T,nbins = 4,full.moderate = F,
                                  ylab = "Marginal Effect of \'Mean Panel Liberalism\'",
                                  xlab = "Moderator: ln(Workload)",
                                  vcov.type = "cluster",
                                  theme.bw = T)
  
  jpeg("../workload_me.jpeg",width=8,height=5, units = "in", res = 800)
  plot(workload_binning,
       theme.bw=T,bin.labs=F,
       xlab="ln(District Court Workload)",
       ylab="Marginal Effect of \'Mean Panel Liberalism\'")
  dev.off()

# table 
  
  stargazer(workload_1,workload_2,workload_3,workload_4,workload_5,
            style = "AJPS",
            keep.stat = c("n", "adj.rsq"),
            header = FALSE,
            covariate.labels = c("Mean Liberalism of Panels",
                                 "ln(Workload)",
                                 "Mean Liberalism $\\times$ Workload"),
            digits = 2,
            label = "workload",
            title = "Predicting liberal district court voting with the mean liberalism of appellate panels, interacted with a measure of district court judges' workloads. Coefficients are from linear fixed effect models with standard errors clustered by circuit.",
            model.names = FALSE,
            add.lines = list(c("Judge Fixed Effects", "", "\\checkmark", "\\checkmark", "\\checkmark", ""),
                             c("Judge-Category Fixed Effects", "", "", "", "", "\\checkmark"),
                             c("Year Fixed Effects", "", "", "\\checkmark", "\\checkmark", "\\checkmark"),
                             c("Case Type Fixed Effects", "", "", "", "\\checkmark", "\\checkmark")),
            dep.var.labels = "Liberal Vote") %>%
    cat(sep = "\n", file = "../workload.tex")

########################
# Progressive Ambition #
########################

promotion_binning <- interflex(data=as.data.frame(district),
                               estimator = "binning",
                               Y="libcon",D="mean_panel",X="time_on_bench",
                                  FE=c("judge_fe","year","casetype"),cl="circuit",
                                  na.rm=T,nbins = 4,full.moderate = F,
                                  ylab = "Marginal Effect of \'Mean Panel Liberalism\'",
                                  xlab = "Time on Bench (Mean-Centered)",
                                  vcov.type = "cluster",
                                  theme.bw = T)

jpeg("../time_me.jpeg",width=8,height=5, units = "in", res = 800)
plot(promotion_binning,
     theme.bw=T,bin.labs=F,
     xlab="District Court Judge Time on Bench (Mean-Centered)",
     ylab="Marginal Effect of \'Mean Panel Liberalism\'")
dev.off()

panel_1_time <- felm(libcon ~ mean_panel * time_on_bench  | 0 | 0 | circuit,
                     data = district,
                     exactDOF = TRUE)

panel_2_time <- felm(libcon ~ mean_panel * time_on_bench  | judge_fe | 0 | circuit,
                     data = district,
                     exactDOF = TRUE)

panel_3_time <- felm(libcon ~ mean_panel * time_on_bench - time_on_bench | judge_fe + factor(year)  | 0 | circuit,
                     data = district,
                     exactDOF = TRUE)

panel_4_time <- felm(libcon ~ mean_panel * time_on_bench - time_on_bench | judge_fe + factor(year) + factor(casetype) | 0 | circuit,
                   data = district,
                   exactDOF = TRUE)

panel_5_time <- felm(libcon ~ mean_panel * time_on_bench - time_on_bench | judge_category + factor(year) + factor(casetype) | 0 | circuit,
                     data = district,
                     exactDOF = TRUE)

panel_1_party <- felm(libcon ~ mean_panel * same_party  | 0 | 0 | circuit,
                     data = district,
                     exactDOF = TRUE)

panel_2_party <- felm(libcon ~ mean_panel * same_party  | judge_fe | 0 | circuit,
                     data = district,
                     exactDOF = TRUE)

panel_3_party <- felm(libcon ~ mean_panel * same_party  | judge_fe + factor(year)  | 0 | circuit,
                     data = district,
                     exactDOF = TRUE)

panel_4_party <- felm(libcon ~ mean_panel * same_party  | judge_fe + factor(year) + factor(casetype) | 0 | circuit,
                     data = district,
                     exactDOF = TRUE)

panel_5_party <- felm(libcon ~ mean_panel * same_party  | judge_category + factor(year) + factor(casetype) | 0 | circuit,
                     data = district,
                     exactDOF = TRUE)

stargazer(panel_1_time,panel_2_time,panel_3_time,panel_4_time,panel_5_time,
          style = "AJPS",
          keep.stat = c("n", "adj.rsq"),
          header = FALSE,
          covariate.labels = c("Mean Liberalism of Panels",
                               "Time on Bench (mean-centered)",
                               "Mean Liberalism $\\times$ Time on Bench"),
          digits = 2,
          label = "time_tab",
          title = "Predicting liberal district court voting with the mean liberalism of appellate panels, interacted with time on bench. Coefficients are from linear fixed effect models with standard errors clustered by circuit.",
          model.names = FALSE,
          add.lines = list(c("Judge Fixed Effects", "","\\checkmark","\\checkmark", "\\checkmark", ""),
                           c("Judge-Category Fixed Effects","","","","","\\checkmark"),
                           c("Year Fixed Effects", "","","\\checkmark", "\\checkmark", "\\checkmark"),
                           c("Case Type Fixed Effects", "","","","\\checkmark", "\\checkmark")),
          dep.var.labels = "Liberal Vote") %>%
  cat(sep = "\n", file = "../time.tex")

stargazer(panel_1_party,panel_2_party,panel_3_party,panel_4_party,panel_5_party,
          style = "AJPS",
          keep.stat = c("n", "adj.rsq"),
          header = FALSE,
          covariate.labels = c("Mean Liberalism of Panels",
                               "Same Party as President",
                               "Mean Liberalism $\\times$ Same Party"),
          digits = 2,
          label = "party_tab",
          title = "Predicting liberal district court voting with the mean liberalism of appellate panels, interacted with co-partisanship with president. Coefficients are from linear fixed effect models with standard errors clustered by circuit.",
          model.names = FALSE,
          add.lines = list(c("Judge Fixed Effects", "","\\checkmark","\\checkmark", "\\checkmark", ""),
                           c("Judge-Category Fixed Effects","","","","","\\checkmark"),
                           c("Year Fixed Effects", "","","\\checkmark", "\\checkmark", "\\checkmark"),
                           c("Case Type Fixed Effects", "","","","\\checkmark", "\\checkmark")),
          dep.var.labels = "Liberal Vote") %>%
  cat(sep = "\n", file = "../party.tex")

# plot for co-partisan with president

party_dat_coefs <- c(panel_4_party$coefficients[1,1],
                   panel_4_party$coefficients[1,1]+panel_4_party$coefficients[3,1])

party_dat_se    <- c(sqrt(panel_4_party$clustervcv[1,1]),
                   sqrt(panel_4_party$clustervcv[1,1]+panel_4_party$clustervcv[3,3]+2*panel_4_party$clustervcv[3,1]))

party_dat <- as.data.frame(cbind(party_dat_coefs,party_dat_se))

party_dat$categories <- c("Not Same Party","Same Party")
party_dat$categories <- factor(party_dat$categories,levels=c("Not Same Party","Same Party"))

jpeg("../party_me.jpeg",width=8,height=5, units = "in", res = 800)
print(
  ggplot(data=party_dat,aes(x=categories,y=party_dat_coefs,ymin=party_dat_coefs-2*party_dat_se,ymax=party_dat_coefs+2*party_dat_se))+
    geom_point(size=4)+
    geom_linerange(size=2)+
    geom_hline(size=2,colour="black",linetype=2,yintercept=0)+
    ylab("Marginal Effect of Mean Panel Liberalism")+
    xlab("District Court Judge is Copartisan with President?")+
    theme_minimal()
)
dev.off()

################################
# logit results ################
################################

logit_1 <- feglm(libcon ~ mean_panel, family = "binomial", data = district)
log1_se <- summary(logit_1,cluster=district$circuit,se="cluster")$se

# Let's try within judge

logit_2 <- feglm(libcon ~ mean_panel  | judge_fe, family = "binomial", data = district)
log2_se <- summary(logit_2,cluster=district$circuit,se="cluster")$se

# Adding year fixed effects

logit_3 <- feglm(libcon ~ mean_panel  | judge_fe + year, family = "binomial", data = district)
log3_se <- summary(logit_3,cluster=district$circuit,se="cluster")$se

# Adding casetype fixed effects

logit_4 <- feglm(libcon ~ mean_panel  | judge_fe + year + casetype, family = "binomial", data = district)
log4_se <- summary(logit_4,cluster=district$circuit,se="cluster")$se

lib_preds <- predict(logit_4,
                     newdata = district %>% mutate(mean_panel = 1))
con_preds <- predict(logit_4,
                     newdata = district %>% mutate(mean_panel = 0))

(lib_preds - con_preds) %>% mean(na.rm = T)

# Judge-category fixed effects

logit_5 <- feglm(libcon ~ mean_panel  | judge_category + year + casetype, family = "binomial", data = district)
log5_se <- summary(logit_5,cluster=district$circuit,se="cluster")$se

logit_placeholder <- lm(libcon ~ mean_panel ,data=district)

stargazer(logit_placeholder,
          logit_placeholder,
          logit_placeholder,
          logit_placeholder,
          logit_placeholder,
          coef = list(coef(logit_1),coef(logit_2),coef(logit_3),coef(logit_4),coef(logit_5)),
          se   = list(log1_se,log2_se,log3_se,log4_se,log5_se),
          style = "AJPS",
          omit.stat = "all",
          header = FALSE,
          covariate.labels = c("Mean Liberalism of Panels"),
          digits = 2,
          label = "logit",
          title = "Predicting liberal district court voting with the mean liberalism of appellate panels, coefficients from logit models (standard errors clustered by circuit)",
          model.names = FALSE,
          add.lines = list(c("Judge Indicators","", "\\checkmark", "\\checkmark", "\\checkmark", ""),
                           c("Judge-Category Indicators", "", "", "", "", "\\checkmark"),
                           c("Year Indicators", "","", "\\checkmark", "\\checkmark", "\\checkmark"),
                           c("Case Type Indicators", "", "","","\\checkmark", "\\checkmark"),
                           c("N",logit_1$nobs,logit_2$nobs,logit_3$nobs,logit_4$nobs,logit_5$nobs),
                           c("Pseudo Adj. $R^2$",round(logit_1$pseudo_r2,3),round(logit_2$pseudo_r2,3)
                             ,round(logit_3$pseudo_r2,3),round(logit_4$pseudo_r2,3),round(logit_5$pseudo_r2,3))),
          dep.var.labels = "Liberal Vote") %>%
  cat(sep = "\n", file = "../logit.tex")


# Correlation between circuit liberalism (median measure) and district judge voting

panel_4_med <- felm(libcon ~ median_all | judge_fe + factor(year) + factor(casetype) | 0 | circuit,
                data = district,
                exactDOF = TRUE)

panel_4_mean <- felm(libcon ~ mean_all| judge_fe + factor(year) + factor(casetype) | 0 | circuit,
                    data = district,
                    exactDOF = TRUE)

stargazer(panel_4_med, panel_4_mean,
          style = "AJPS",
          keep.stat = c("n", "adj.rsq"),
          header = FALSE,
          covariate.labels = c("Median Liberalism of Circuit Judges",
                               "Mean Liberalism of Circuit Judges"),
          label = "felm_median",
          title = "Predicting liberal district court voting with alternative measures of circuit court liberalism; coefficients from linear fixed effect models (standard errors clustered by circuit)",
          digits = 2,
          model.names = FALSE,
          add.lines = list(c("Judge Fixed Effects","\\checkmark","\\checkmark"),
                           c("Year Fixed Effects","\\checkmark","\\checkmark"),
                           c("Case Type Fixed Effects","\\checkmark","\\checkmark")),
          dep.var.labels = "Liberal Vote") %>%
  cat(sep = "\n", file = "../felm_median.tex")

# alternative levels of measurement
jy_data <- district %>% group_by(judge_fe,statdist,circuit,year) %>% summarise(libcon=mean(libcon,na.rm=T),mean_panel=mean(mean_panel,na.rm=T))
dy_data <- district %>% group_by(statdist,circuit,year) %>% summarise(libcon2=mean(libcon,na.rm=T),mean_panel=mean(mean_panel,na.rm=T))
cy_data <- district %>% group_by(circuit,year) %>% summarise(libcon3=mean(libcon,na.rm=T),mean_panel=mean(mean_panel,na.rm=T))
  
panel_3_jy <- felm(libcon ~ mean_panel | judge_fe + factor(year) | 0 | circuit, data = jy_data)
panel_3_dy <- felm(libcon2 ~ mean_panel | statdist + factor(year) | 0 | circuit, data = dy_data)
panel_3_cy <- felm(libcon3 ~ mean_panel | circuit + factor(year) | 0 | circuit, data = cy_data)




stargazer(panel_3_jy, panel_3_dy, panel_3_cy,
          style = "AJPS",
          keep.stat = c("n", "adj.rsq"),
          header = FALSE,
          covariate.labels = c("Mean Liberalism of Circuit Panels"),
          label = "felm_levels",
          title = "Predicting liberal district court voting with the mean liberalism of appellate panels; coefficients from linear fixed effect models at different levels of aggregation (standard errors clustered by circuit)",
          digits = 2,
          model.names = FALSE,
          add.lines = list(c("Judge Fixed Effects", "\\checkmark", "", ""),
                           c("District Fixed Effects", "", "\\checkmark", ""),
                           c("Circuit Fixed Effects", "", "", "\\checkmark"),
                           c("Year Fixed Effects", "\\checkmark", "\\checkmark", "\\checkmark")),
          dep.var.caption = "Liberal Vote",
          dep.var.labels = c("Judge-Year Level","District-Year Level","Circuit-Year Level")) %>%
  cat(sep = "\n", file = "../felm_levels.tex")

# scatter-plot at circuit-year level

  cor_out <- cor(cy_data$libcon3,cy_data$mean_panel,use="complete.obs")

  fileConn<-file("../cor_out.tex")
  writeLines(format(round(cor_out,3)), fileConn)
  close(fileConn)

  jpeg("../circ_level_scatter.jpeg",width=6,height=5, units = "in", res = 800)
  print(
    ggplot(data=cy_data,aes(y=libcon3,x=mean_panel))+
      geom_point(colour="grey20",alpha=0.4)+
      geom_smooth(size=2,se=F,colour="black",method="lm")+
      theme_minimal()+
      xlab("Circuit Court Mean Panel Liberalism")+
      ylab("Proportion District Court Decisions in Liberal Direction")
  )
  dev.off()
 

# alternative specifications
 
  district$year_counter <- as.numeric(district$year)-1964
  cy_data$year_counter  <- as.numeric(cy_data$year)-1964
  cy_data <- dplyr::arrange(cy_data,circuit,year)
  cy_data <- slide(cy_data,Var="libcon3",NewVar="ldv",GroupVar="circuit",TimeVar="year",slideBy = -1)
  
  panel_3_lt <- felm(libcon ~ mean_panel | as.numeric(year_counter):factor(circuit) + factor(judge_fe) + factor(year) + factor(casetype)| 0 | circuit, data = district) 

  panel_3_cy_lt <- felm(libcon3 ~ mean_panel | as.numeric(year_counter):factor(circuit) + factor(circuit) + factor(year)| 0 | circuit, data = cy_data) 

  panel_3_cy_ldv <- felm(libcon3 ~ mean_panel + ldv| year | 0 | circuit, data = cy_data) 
  
  panel_3_cy_ldv_fe <- felm(libcon3 ~ mean_panel + ldv| circuit + year | 0 | circuit, data = cy_data) 
  
  stargazer(panel_3_lt,panel_3_cy_lt,panel_3_cy_ldv,panel_3_cy_ldv_fe,
            style = "AJPS",
            keep.stat = c("n", "adj.rsq"),
            header = FALSE,
            covariate.labels = c("Mean Liberalism of Circuit Panels","Lagged DV"),
            label = "altspecs",
            title = "Predicting liberal district court voting with the mean liberalism of appellate panels; coefficients from linear models with alternative methods of controlling for unobservables (standard errors clustered by circuit)",
            digits = 2,
            model.names = FALSE,
            add.lines = list(c("Judge Fixed Effects", "\\checkmark", "","",""),
                             c("Circuit Fixed Effects", "", "\\checkmark", "", "\\checkmark"),
                             c("Circuit Linear Trends", "\\checkmark", "\\checkmark", "",""),
                             c("Year Fixed Effects", "\\checkmark", "\\checkmark", "\\checkmark", "\\checkmark")),
            dep.var.labels = c("Judge-Decision","Circuit-Year")) %>%
    cat(sep = "\n", file = "../felm_altspecs.tex")

###################################################################################################
# Leave-One-Out Analyses ##########################################################################
###################################################################################################
 
# circuits 
 
 circ_list <- list()
 out_circuits <- list()
 
 for(j in 1:length(unique(district$circuit))){
   
   circuit <- sort(unique(district$circuit))[j]
   
   dat_circ <- district[district$circuit!=circuit,]
   
   panel_4_circ <- felm(libcon ~ mean_panel | judge_fe + factor(year) + factor(casetype) | 0 | circuit,
                   data = dat_circ,
                   exactDOF = TRUE)
   
   circ_list[[j]] <- summary(panel_4_circ)$coefficients[1,]
   out_circuits[[j]]   <- circuit
   
 }
 
 circ_plot_dat <- as.data.frame(do.call(rbind,circ_list)); rownames(circ_plot_dat) <- NULL
 circ_plot_dat$Circuit <- unlist(out_circuits)
 circ_plot_dat$Circuit <- factor(circ_plot_dat$Circuit,levels=c("First","Second","Third","Fourth",
                                                                "Fifth, Old","Fifth, New","Sixth",
                                                                "Seventh","Eighth","Ninth","Tenth",
                                                                "Eleventh","D.C."))
 
 jpeg("../circuit_out.jpeg",width=10,height=5, units = "in", res = 800)
 print(
   ggplot(data=circ_plot_dat,aes(x=Circuit,y=Estimate,ymin=Estimate-2*`Cluster s.e.`, ymax=Estimate+2*`Cluster s.e.`))+
     geom_hline(yintercept=0,size=2,colour="indianred",linetype=2)+
     geom_point(size=4)+
     geom_linerange(size=2)+
     theme_minimal()+
     ylab("Estimate of 'Mean Panel Liberalism'")
 )
 dev.off()
 
# year
 
 district$decade <- str_sub(district$year,1,3)
 
 dec_list <- list()
 out_dec <- list()
 
 for(j in 1:length(unique(district$decade))){
   
   decade <- sort(unique(district$decade))[j]
   
   dat_dec <- district[district$decade!=decade,]
   
   panel_4_dec <- felm(libcon ~ mean_panel | judge_fe + factor(year) + factor(casetype) | 0 | circuit,
                        data = dat_dec,
                        exactDOF = TRUE)
   
   dec_list[[j]] <- summary(panel_4_dec)$coefficients[1,]
   out_dec[[j]] <- decade
   
 }
 
 dec_plot_dat <- as.data.frame(do.call(rbind,dec_list)); rownames(dec_plot_dat) <- NULL
 dec_plot_dat$Decade <- paste(unlist(out_dec),"0s",sep="")
 dec_plot_dat$Decade <- factor(dec_plot_dat$Decade,levels=paste(sort(unique(district$decade)),"0s",sep=""))

 jpeg("../decade_out.jpeg",width=8,height=5, units = "in", res = 800)
 print(
   ggplot(data=dec_plot_dat,aes(x=Decade,y=Estimate,ymin=Estimate-2*`Cluster s.e.`, ymax=Estimate+2*`Cluster s.e.`))+
     geom_hline(yintercept=0,size=2,colour="indianred",linetype=2)+
     geom_point(size=4)+
     geom_linerange(size=2)+
     theme_minimal()+
     ylab("Estimate of 'Mean Panel Liberalism'")
 )
 dev.off()
 
# issue areas 
 
 district$issue <- NA
 district$issue <- as.character(district$issue)
 district[district$casetype==1,"issue"] <- "Habeus Corpus - U.S."
   district[district$casetype==2,"issue"] <- "Habeus Corpus - State"
   district[district$casetype==3,"issue"] <- "Criminal Court Motion"
   district[district$casetype==4,"issue"] <- "Contempt of Court"
   district[district$casetype==5,"issue"] <- "Non-Conv. Criminal"
   district[district$casetype==6,"issue"] <- "Alien Petitions"
   district[district$casetype==7,"issue"] <- "Native American Rights"
   district[district$casetype==8,"issue"] <- "Voting Rights"
   district[district$casetype==9,"issue"] <- "Social Security"
   district[district$casetype==10,"issue"] <- "Racial Discrimination"
   district[district$casetype==11,"issue"] <- "14th Amendment"
   district[district$casetype==12,"issue"] <- "Military Exclusion"
   district[district$casetype==13,"issue"] <- "Freedom of Expression"
   district[district$casetype==14,"issue"] <- "Freedom of Religion"
   district[district$casetype==15,"issue"] <- "Union vs. Company"
   district[district$casetype==16,"issue"] <- "Member vs. Union"
   district[district$casetype==17,"issue"] <- "Employee vs. Employer"
   district[district$casetype==18,"issue"] <- "Commercial Regulation"
   district[district$casetype==19,"issue"] <- "Environmental Protection"
   district[district$casetype==20,"issue"] <- "Local/State Economic"
   district[district$casetype==21,"issue"] <- "Govt. vs. Union/Employer"
   district[district$casetype==22,"issue"] <- "Rent Control/Excess Profits"
   district[district$casetype==23,"issue"] <- "Women's/Gender Rights"
   district[district$casetype==24,"issue"] <- "NLRB vs. Employer"
   district[district$casetype==25,"issue"] <- "NLRB vs. Union"
   district[district$casetype==26,"issue"] <- "Disabled Rights"
   district[district$casetype==27,"issue"] <- "Reverse Discrimination - Race"
   district[district$casetype==28,"issue"] <- "Reverse Discrimination - Sex"
   district[district$casetype==29,"issue"] <- "Right to Privacy"
   district[district$casetype==30,"issue"] <- "Age Discrimination"
   district[district$casetype==31,"issue"] <- "Sentencing Guidelines Deviation"
     
 iss_list <- list()
 out_issue <- list()
 
 for(j in 1:length(unique(district$issue))){
   
   issue <- sort(unique(district$issue))[j]
   
   dat_iss <- district[district$issue!=issue,]
   
   panel_4_iss <- felm(libcon ~ mean_panel | judge_fe + factor(year) + factor(casetype) | 0 | circuit,
                        data = dat_iss,
                        exactDOF = TRUE)
   
   iss_list[[j]] <- summary(panel_4_iss)$coefficients[1,]
   out_issue[[j]] <- issue
   
 }
 
 iss_plot_dat <- as.data.frame(do.call(rbind,iss_list)); rownames(iss_plot_dat) <- NULL
 iss_plot_dat$Issue <- unlist(out_issue)
 iss_plot_dat$Issue <- factor(iss_plot_dat$Issue,levels=sort(unique(district$issue)))
 
 jpeg("../issue_out.jpeg",width=8,height=5, units = "in", res = 800)
 print(
   ggplot(data=iss_plot_dat,aes(x=Issue,y=Estimate,ymin=Estimate-2*`Cluster s.e.`, ymax=Estimate+2*`Cluster s.e.`))+
     geom_hline(yintercept=0,size=2,colour="indianred",linetype=2)+
     geom_point(size=4)+
     geom_linerange(size=2)+
     theme_minimal()+
     theme(axis.text.x=element_text(angle = 45, vjust=1,hjust=1))+
     ylab("Estimate of 'Mean Panel Liberalism'")
 )
 dev.off()
 

 ## Cohort Time Trends
 
 
cohort_time <- felm(libcon ~ mean_panel | factor(judge_fe) + factor(appres):as.numeric(year-2000) + factor(year) + factor(casetype) | 0 | circuit,
      data = district,
      exactDOF = TRUE)

stargazer(cohort_time,
          style = "AJPS",
          keep.stat = c("n", "adj.rsq"),
          header = FALSE,
          covariate.labels = c("Mean Liberalism of Circuit Panels"),
          label = "cohort_time",
          title = "Predicting liberal district court voting with the mean liberalism of appellate panels; coefficients from linear fixed effect models including time trends by appointing president (standard errors clustered by circuit)",
          digits = 2,
          model.names = FALSE,
          add.lines = list(c("Judge Fixed Effects", "\\checkmark"),
                           c("Year Fixed Effects", "\\checkmark"),
                           c("Case Type Fixed Effects", "\\checkmark"),
                           c("Appointing President Time Trends", "\\checkmark")),
          dep.var.labels = "District Court Judge Liberal Vote") %>%
  cat(sep = "\n", file = "../cohort_time.tex")
